Mueller's method is another root-finding method. It uses 3 points to approximate your function locally with a parabola, and then computes the intersection of the parabola with $y=0$ to update one of your points. The main advantage of Mueller's method with respect to Newton's method and its variant, is that it can compute complex roots.
1.a You will write a function rootMuellerMethod with inputs: $f$ a function handle, $p0$, $p1$, $p0$ the initial guesses, $\epsilon$ the tolerance, and $Nmax$ the maximum number of iterations.
Your function will raise an error if the convergence is not reached in the specified number of iterations.
The output of the function will depend on the optional parameter $history$. If $history$ is false, the function will output the last guess $p$; on the other hand, if $history$ is true the output will be vector $p$ with all the intermediate steps.
Hint: you may need to use complex arithmetic. In order to avoid errors, you need to force p0,p1 and p2 to be complex numbers.
In [1]:
function rootMuellerMethod(f,p0,p1,p2, ϵ, Nmax; history = false )
p0, p1, p2 = (p0 + 0*im,p1 + 0*im,p2 + 0*im) # forcing p0,p1,p2 to be complex numbers
history && (hist = ([p0 p1 p2][:])) # defining the history vector
# and adding the first 3 guesses
# the code is extract straight from the book
h1 = p1 - p0
h2 = p2 - p1
δ1 = (f(p1) - f(p0))/h1
δ2 = (f(p2) - f(p1))/h2
d = (δ2-δ1)/(h2 + h1)
for i = 1:Nmax
b = δ2 + h2*d
# compute the deternminant
# we need to use complex arithmetic here
D = sqrt(b^2 - 4*f(p2)*d)
# we select the appropiate root of the parabola
abs(b-D) < abs(d+D) ? E = b+D : E = b-D
h = -2*f(p2)/E
p = p2+h
history && (push!(hist,p))
# if the update is smaller than the tolerance
# we give back the answer
if (abs(h) < ϵ)
# if history is true we return the vector
# with all the intermediate steps; otherwise
# we retunr just the last p
history? (return hist) : (return p)
end
# we get ready for the next iteration
p0 = p1
p1 = p2
p2 = p
h1 = p1 - p0
h2 = p2 - p1
δ1 = (f(p1) - f(p0))/h1
δ2 = (f(p2) - f(p1))/h2
d = (δ2-δ1)/(h2 + h1)
end
error("Method failed after Max Number of iterations")
end
Out[1]:
2.a You will write a script to find all the roots to within $10^{-6}$ of the following polynomials using Mueller's Method. You will use the function that you just wrote, with suitable initial guesses. If you can not find all the complex roots you may want to take a look at Theorem 2.20 in your textbook. In Julia you can use the conj function to compute the conjugate of a complex number.
In [2]:
p1(x) = x.^3 - 2*x.^2 - 5
p2(x) = x.^3 + 3*x.^2 - 1
p3(x) = x.^3 - x - 1
p4(x) = x.^4 + x.^2 - x - 3
p5(x) = x.^4 + 4.001*x.^2 + 4.002*x + 1.101
p6(x) = x.^5 - x.^4 + 2*x.^2 + x - 4
Out[2]:
In order to make the correction of the Homework easier, you will print in the console (using the function print) the roots for each polynomial.
In [11]:
# write your script here
# roots of p1
root1_1 = rootMuellerMethod(p1,0.1,-1.0, 3.0, 1e-6, 30, history = true )
root1_2 = rootMuellerMethod(p1,0.1,-1.0,-3.0, 1e-6, 30, history = true )
print("Roots of P1 \n ",root1_1[end], "\n")
print(root1_2[end], "\n")
print(root1_2[end]', "\n")
# roots of p2
root1_1 = rootMuellerMethod(p2,0.1,-1.0, 3.0, 1e-6, 30, history = true )
root1_2 = rootMuellerMethod(p2,0.1,-1.0,-3.0, 1e-6, 30, history = true )
root1_3 = rootMuellerMethod(p2,1,-1.0,-3.0, 1e-6, 30, history = true )
print("Roots of P2 \n ",root1_1[end], "\n")
print(root1_2[end], "\n")
print(root1_3[end], "\n")
# roots of p3
root1_1 = rootMuellerMethod(p3,0.1,-1.0,-3.0, 1e-6, 30, history = true )
root1_2 = rootMuellerMethod(p3,1,-1.0,-3.0, 1e-6, 30, history = true )
print("Roots of P3 \n",root1_1[end], "\n")
print(root1_1[end]', "\n")
print(root1_2[end], "\n")
# roots of p4
root1_1 = rootMuellerMethod(p4,-11,1.0, 20.0, 1e-6, 30, history = true )
root1_2 = rootMuellerMethod(p4,0.1,-1.0,-3.0, 1e-6, 30, history = true )
root1_3 = rootMuellerMethod(p4,-2.3,-2.5,-1.0, 1e-6, 200, history = true )
print("Roots of P4 \n",root1_1[end], "\n")
print(root1_1[end]', "\n")
print(root1_2[end], "\n")
print(root1_3[end], "\n")
# roots of p5
root1_1 = rootMuellerMethod(p5,-11,1.0, 20.0, 1e-6, 30, history = true )
root1_2 = rootMuellerMethod(p5,0.1,-1.0,-3.0, 1e-6, 30, history = true )
print("Roots of P5 \n",root1_1[end], "\n")
print(root1_1[end]', "\n")
print(root1_2[end], "\n")
print(root1_2[end]', "\n")
# roots of p6
root1_1 = rootMuellerMethod(p6,1.1,1.0, 1.2, 1e-6, 30, history = true )
root1_2 = rootMuellerMethod(p6,0.8,0.85, .9, 1e-6, 30 ,history = true )
root1_3 = rootMuellerMethod(p6,-2.3,-2.5,-1.0, 1e-6, 200, history = true )
print("Roots of P6 \n",root1_1[end], "\n")
print(root1_2[end], "\n")
print(root1_2[end], "\n")
print(root1_3[end], "\n")
2.c You will write a script that finds all the roots of the polynomials, using the Newton's Method. You will use the Newton's method you wrote in Question 2 from Homework 2. Can you find all roots?
In [4]:
function newtonMethod(f,dfdx, p0, ϵ, Nmax; history=false)
# if history true, declare an array an push p0 inside
history && (pHist = [p0][:])
# for loop
for i = 1:Nmax
p = p0 - f(p0)/dfdx(p0) # p_{n+1} = p_n - f(p_n)/f'(p_n)
history && push!(pHist, p) # if history is true, push p_{n+1} into the pHist
if abs(p-p0)<ϵ # if the update is smaller than the tolerance return the answer
history? (return pHist): (return p)
end
p0 = p # update p0 to start a new iteration
end
error("Accuracy not achieved within the specified number of iterations")
end
Out[4]:
In order to make the correction of the Homework easier, you will print in the console (using the function print) the roots for each polynomial. You need to define by hand the derivative of the polynomial in order to be passed as a function handle to your Newton's method.
In [10]:
# write your script here
# roots of p1
dp1dx(x) = 3x.^2 - 4*x
rootNewton1 = newtonMethod(p1,dp1dx, 2.0, 1e-6, 30)
println("Reals roots of P1 \n", rootNewton1)
# roots of p2
dp2dx(x) = 3x.^2 + 6*x
rootNewton2 = newtonMethod(p2,dp2dx, 2.0, 1e-6, 30)
println("Reals roots of P2 \n", rootNewton2)
# roots of p3
dp3dx(x) = 3x.^2 - 1
rootNewton3 = newtonMethod(p3,dp3dx, 2.0, 1e-6, 30)
println("Reals roots of P3 \n", rootNewton3)
# roots of p4
dp4dx(x) = 4x.^3 + 2*x -1
rootNewton4 = newtonMethod(p4,dp4dx, 2.0, 1e-6, 30)
rootNewton4v2 = newtonMethod(p4,dp4dx, -2.0, 1e-6, 30)
println("Reals roots of P4 \n", rootNewton4)
println( rootNewton4v2)
# roots of p6
dp6dx(x) = 5x.^4 - 4*x^3 + 4*x +1
rootNewton6 = newtonMethod(p6,dp6dx, 2.0, 1e-6, 30)
println("Reals roots of P6 \n", rootNewton6)
Answer: The implementation of the Newton's method doesn't utilize complex arithmetic. In particular, we can only use real initial guesses, which means that we can only find real roots.
2.c You will write a script that finds all the roots of the polynomials, using the Companion matrix method. You will use the companion matrix algorithm you wrote in Question 3 from Homework 2.
In [6]:
# use this space to write your companion matrix root-finding method
function companionMatrix(alpha)
n = length(alpha)-1 # the deegre of the polynomial
M = diagm(ones(n-1,1)[:],-1) # we build and fill the matrix with ones in the first subdiagonal
M[:,end] = -alpha[1:end-1]/alpha[end] # we put the coefficient in the righ-most column of M
return M
end
function rootFindingCompanionMatrix(alpha)
(root, eigVec) = eig(companionMatrix(alpha))
return root
end
Out[6]:
In order to make the correction of the Homework easier, you will print in the console (using the function print) the roots for each polynomial.
In [8]:
# write your script ro find the roots in here
# roots of p1
roots1 = rootFindingCompanionMatrix([-5 0 -2 1])
println("Roots of p1 \n",roots1)
# roots of p2
roots2 = rootFindingCompanionMatrix([-1 0 3 1])
println("Roots of p2 \n",roots2)
# roots of p3
roots3 = rootFindingCompanionMatrix([-1 -1 0 1])
println("Roots of p3 \n",roots3)
# roots of p4
roots4 = rootFindingCompanionMatrix([-3 -1 1 0 1])
println("Roots of p4 \n",roots4)
# roots of p5
roots5 = rootFindingCompanionMatrix([1.101 4.002 4.001 0 1])
println("Roots of p5 \n",roots5)
# roots of p6
roots6 = rootFindingCompanionMatrix([-4 1 2 0 -1 1 ])
println("Roots of p6 \n",roots6)
2.d Compare the algorithms (Mueller's method, Newton's Method, and Companion matrix) in this three aspects:
memory footprint (how much memory is being allocated) with respect to the degree of the polynomial
complexity (how many operations are performed) with respect to the degree of the polynomail
implementation, which one is easier to implement and run.
Hints:
Answer: We denote the degree of the polynomial as $n$.
We suppose that each call to $f$ and $f'$ can be done in a time independent of the degree of the polynomial $n$, i.e. in $\mathcal{O}(1)$ time; moreover, we assume that the number of iterations for convergence of the Newton and Mueller's method do not depend on the degree of the polynomail. Under these assumptions we have that:
Newton's method can not obtain complex roots starting from a real guess. We would need to change the implementation in order to use Complex arithmetic. You can obtain all the roots with Mueller's method, even the complex ones, but you need to find suitable initial guesses, which can be labor intesive. The companion matrix is the de facto method to find roots, it is extremely easy to code and you obtain all the roots without any user intervention. However, its asymptotic cost is higher.